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Abstract 


We present a new theory for modeling forced indentation spectral lineshapes of biological parti¬ 
cles, which considers non-linear Hertzian deformation due to an indenter-particle physical contact 
and bending deformations of curved beams modeling the particle structure. The bending of beams 
beyond the critical point triggers the particle dynamic transition to the collapsed state, an extreme 
event leading to the catastrophic force drop as observed in the force (-F)-deformation (A) spectra. 
The theory interprets fine features of the spectra: the slope of the FX curves and the position of 
force-peak signal, in terms of mechanical characteristics — the Young’s moduli for Hertzian and 
bending deformations En and £),, and the probability distribution of the maximum strength with 
the strength of the strongest beam and the beams’ failure rate rn. The theory is applied to 
successfully characterize the FX curves for spherical virus particles — CCMV, TrV, and AdV. 
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Mechanical testing has become the main tool to probe the physico-chemical and mate 


Q. 


rials properties of the protein shells of plant and animal viruses, and bacteriophages 
A variety of viruses have been explored by profiling the indentation force F as a function 
of particle deformation X (.FX curve), including bacteriophages $29 and HK97 ffl, the 
human viruses Noro Virus, Hepatitis B Virus, Human Immuno Deficiency Virus (HIV), Ade¬ 


novirus 


(AdV) and Herpes Simplex Virus [5j-8] and other eukaryotic cell infecting 


viruses 


like Minute Virus of Mice, 1 
Virus (CCMV) and BMV 


’ riatoma Virus (TrV) and plant viruses Cowpea Chlorotic Mottle 
9-113|. Although these experiments reveal an amazing diversity 
of mechanical properties of biological particles, experimental results are difficult to inter¬ 
pret without theoretical modeling. What types of mechanical excitations drive the particle 
deformation and collapse? What determines the mechanical limit of the particle — the crit¬ 
ical forces and critical deformations? Why is the initial portion of the FX spectra weakly 
non-linear? Why do the FX spectra for the same particle differ from one measurement to 
another? This points to the stochastic nature of collapse transitions, but what defines the 
likelihood of structural collapse at a given force load? 

A number of theoretical approaches have been d esig ned to describe the dynamics of 
virus particles, including: finite element analysis jl4, 15], normal mode analysis 
tic network modeling 17J, atomistic MD and coarse-grained simulations 18 


16], elas- 


approaches 


20], and other 


21 ], Here we take a step further to develop an analytically tractable model for 


meaningful interpretation of the force-deformation spectral lineshapes available from single- 
particle nanomanipulation experiments. The theory links the slope, critical force, and the 
critical deformation of the FX curve with the physical characteristics of the structure, ge¬ 
ometry and overall shape of the particle and indenter. We identify the types of mechanical 
excitations (degrees of freedom), which contribute to the particle deformation (indentation 
depth) X, by analyzing the structure and potential energy changes in the CCMV particle 
using our results of nanoindentations in silico. The methodology of nanoindentation in silico 
is described in Ref. [22]; see Supplementary Material (SM; Figs. S1-S3). We formulate and 
apply the Fluctuating Nonlinear Spring (FNS) model to characterize the FX spectra for 
the CCMV, AdV and TrV particles obtained as described in Refs.jg, ll|, 23]. 

Degrees of freedom — In dynamic force-ramp assays f(t)—KUft, an indenter (cantilever 
tip) compresses a particle (Fig. Q] S3) creating a physical contact between them. The force 
loads the particle mechanically over time t with the force-loading rate kvj (k and Vf are 
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the cantilever spring constant and velocity). For small /, the mechanical energy is localized 
to the virion surface under the tip, which results in the tip and particle undergoing nor¬ 
mal displacements u tip and u par corresponding to the Hertzian deformation XH=u t i P +u par . 
Typically, u tip <^u par , therefore XH=u par . The mechanical energy gradually fills the parti¬ 
cle structure stressing the side portions of the structure undergoing bending deformations 
Xb (Fig. [[]). The total deformation is the sum of Hertzian and bending deformations, i.e. 
X=XH-\-Xb, and the deformation force F of the particle (restoring force exerted on the tip 
by the particle) is a bivariate function, F (X=x H+Xb)=FH(;x h)+ F b (;Xb) ■ We quantified the 
amplitude of Hertzian deformation Xh ~3 nm and bending deformations x b ~ 4.3 nm using 
the simulation output for the CCMV (Fig. [T]). 


The Hertz model 
(Fig.[U(a)), 


24 


25] accounts for the force due to local curvature change F h (xh ) 


Fh(xh) — 1 /DniRparRtip/iRpar + Etip)Y^ ' F’h 


( 1 ) 


where R par and R Up are the radii of the particle and the tip and D H =0.7b((l — o 2 H ) / E H + (l — 
a ti P ) / E t i P ) ■ E h and E tip are the Young’s moduli and an and a tip are the Poisson’s ratios for 
the particle and the tip deformations, respectively. Since Eup^En, Dh—0-75(1 — o 2 h )/Eh- 
To describe the bending deformations Fb(x b ), we divide the side portion of the particle 
structure into vertical curved beams of length L (Fig. [0(d)). For a spherical particle with 
thickness r, the total number of beams is N=2nR par /r. For small x b (Fig. [T]), the potential 
energy change is E b I /2 f L (K,(xb,l)—Ko) 2 dl 24, [26|, where Ko and K(x bl l ) are the initial and 
current curvature of the beam element dl (0<l<R par — Xb/2) and EbI is the flexural rigidity. 
With the beam shape function q(xb , l) = (R par +Xb/ 2)(1 — l 2 /( Rp ar — Xb/ 2) 2 ) 1 / 2 the curvature 
of the beam is given by n(xb, l)=q"(xb, l)/(l+(q'(xb, l)) 2 ) 3//2 (q' and q" are the first and 
second derivatives of q with respect to l). By performing the integration we obtain the 
expression for the beam bending energy, which upon differentiation with respect to Xb, gives 
the beam-bending force. Expanding the resulting expression in Taylor series in powers of Xb 
and retaining the linear term we obtain: 


fb{xb) = 9E b lTr/8Rl ar ■ x b (2) 

Combining the contributions from all N beams and adding together Eqs.([T]) and (121) , 
we obtain the deformation force F(xh, Xb)=FH(xH)+Fb(xb)=kHX 3 Ip 2 +NkbXb 1 where kn = 
(R par Rti P /(R P ar + Rtip) 1 ^ 2 )/Dh (Eq.flT])) and kb—9E b In/(8R par ) is the beam spring constant. 
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In agreement with the results of in silico indentation of CCMV (Fig. [2]), F(xH,x b ) predicts 
that the initial portion of FX curves, where Hertzian deformation dominates (Xpzxh), is 


weakly nonlinear. The Hertzian nonlinearity expressed in F(xjj,x b ) was indeed observed in 


recent experiments on thick-shelled particles 10]. However, F(xH,x b ) does not capture the 


catastrophic force drop (Fig. [2]) because the theory lacks a description of structural damage. 


Fluctuating Nonlinear Spring (FNS) model — We represent a particle by a collection of 
N identical beams interacting with an indenter through a Hertzian cushion. Each i-th beam 
undergoes the equal elastic deformation Xbi=x b with the spring constant (z=l, 2,..., N) 
until it fails mechanically when the load on the beam reaches some critical value / b *. The 
spherical geometry of a virus particle dictates the parallel arrangement of beams with the 
spring K b =Eti kbi=Nk b . There are n (or N—n) beams that have failed (or survived), and 
the actual bending force is Fb(xb)=kb(N — n)x b —K b x b ( 1 — n/N). We define the random 
variables associated with the probability of damage d—n/N and survival s={N — n)/N=l — d 
of the collection of beams (0<d, s<l). In the continuous limit, d and s are described by the 


pdfs p(F b ), d(F b )=Prob(F b =K b x b )=f 0 Ft p(F£)dF£ and s{F b )= 1 - d{F b )=J%p{FfidFl. With 


the particle damage accounted for, the deformation force becomes: 


F(x H ,x b ) = k H x + K b x b s(x b ) 


(3) 


Here, s can be estimated from computer simulations using the structural similarity between 
a given structure and the initial state x( x b)=(2N(N — l))^ 1 X) ©(| Tij(x b ) — r^(0)| — /3r^(0)). 
In the Heaviside step function, rij(x b ) and r^(0) are the distances between the i-th and j- 
th amino acids in the given and initial structure, respectively, and /3=0.2 is the tolerance 
for the distance change. Because x{x b ) ranges from y—1 (fully similar structures) to y=0 
(completely dissimilar structures) and changes in s(x b ) are commensurate with structure 
changes, we have s(x b )~x(x b ). The y-based estimate of s(x b ) from force-deformation in 
silico of CCMV shows that this metric can be used to inform the modeling of p(F b ) (Fig. [2]). 

The transition to the collapsed state occurs when all beams have failed and so, the longest 
lasting beam determines the collapse onset at the critical force F co1 (critical deformation 
X co1 ). Hence, the statistics of the extreme (maximum) force determines the beams’ failure. 


We used the Weibull distribution 



d(x b ) = 1 - S ( X6 ) 


(4) 
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with the shape parameter m (failure rate), and scale parameter F b which can be under¬ 
stood by using the condition of maximum force. By substituting s(x b ) into Fb(x b ) in 
Eq.pD, by differentiating this expression with respect to x b , and by setting dF b /dx b = 0, 
we obtain Ff —K b x* h y/rn (x b is the critical beam deformation). By substituting Ff into 
the expression for F b (x b ) and into Eq. @ . we obtain the beam-bending force threshold 
pcoi—F*/ v^em=K b x b / \fe and the critical value of collapse probability d col =1 — s col =l — 
1/ ^e. 

By substituting Eq.Q for s(x b ) in Eq.()3|). we obtain the main result of this Letter: 

F(x h , x b ) = k H x 3 H 2 + K b x b e~( KbXb/F b ) " 1 (5) 

Eq.()5ll shows that a biological particle behaves as a nonlinear spring. The beams’ bending 
starts as elastic ( Nk b ), but becomes stochastic near the collapse transition, where the par¬ 
ticle mechanical resistance fluctuates, thus explaining the variability of F co1 and X co1 in the 
FX spectra (Figs. [2j [3]). Hence, the uniaxial forced deformation of a biological particle can 
be represented by the mechanical evolution of a fluctuating nonlinear spring (FNS). 

Application of FNS model — In the experiment, F is measured as a function of the sum 
X=xh + Xb- A particular realization of the forced deformation process (FX trajectory) is 
a stochastic path on the 2D surface F(xH,x b ) (Fig. EJ). For slow loading, when the particle 
structure equilibrates on a timescale faster than the rate of force change, the most dominant 
path is the equilibrium path F eq (X eq ). Using slow cantilever velocities (vf— 0.06—1.0 /irn/s) 
allows us to use this quasi-equilibrium argument. The equilibrium force F eq can then be de¬ 
termined from the requirement that the deformation force (deformation energy) attains the 
minimum. Finding F eq is equivalent to Ending a pair xh and x b that minimizes F(xn,x b ) 
subject to the constraint, X—xh + x b (Fig. [3]), which can be solved using Lagrange multi¬ 
pliers (SM). 

The average simulated spectra for CCMV are compared with the theoretical curves in 
Fig. [2] (simulated force-deformation spectra for the CCMV particle are in Fig. S4(a)). 
To find the best fit, we employed two methods. The exact method is based on Eq. (EJ) 
and uses Lagrange multipliers to find Xh and x b - This requires solving for x b the non¬ 
linear equation aixf m + a 2 a;f m + a 3 xl m + a±x™ + a 5 x b + a 6 = 0, where a\ =rriK^(K b / Ff) 4rn , 
a 2 =-2m(l + m)K 2 b (K b /F b *) 3m , a 3 =(l + 4m + m 2 )K 2 b (K b /F b *) 2m , a 4 =-2(l + m)K 2 {K b /F b *) m , 
a 5 =9/Akjj, and a 6 =Ii 2 — 9/4 k\X. Then, xh is obtained as xn—X—x b . In the piece- 
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wise approximate method (Fig. [2]), the spectrum is divided into the Hertzian-deformation- 
dominated initial regime I: X^xh (xb~0) and F^FH—knX 3 ^ 2 -, and the transition regime 
II (non-monotonic part of FX curve): X^x b and F^F b =K b Xs(X). We calculate F H in 
regime I for XmxH<x^ ax ', Xff ax is obtained using Lagrange multipliers and setting s(a;&)=l. 
In regime II, we use F(Xff ax , x b )=k H (x 7 ff ax ) 3 ^ 2 + K b (x b — x™ x exp [— (K b (x b — x^ ax ) / F£) m }) 
for X^x b >x^ ax . 

The agreement between the simulated force-deformation spectra and theoretical FX 
curves for the CCMV particle is almost quantitative (Fig. [2]); model parameters obtained 
using both methods are very close (Table I). For all symmetry types, the Hertzian excitation 
is softer than the beam-bending ( ku<K b ) implying smaller Young’s modulus for Hertzian 
deformation, Ex<E bl which is why the Hertzian degree of freedom is excited first (Fig. [[]). 
After the Hertzian force reached the maximum F^ ax =kH{x 1 Jf ,x ) 3 ^ 2 at X^XH=x ]} ax , a sub¬ 
sequent force increase excites the beam-bending degrees of freedom and % ( x b ) decreases 
(increases). Hence, physical properties of the particle are dynamic as the nature of its me¬ 
chanical response changes with increasing A" from Hertzian to beam-bending deformation, 
because the actual stiffness of beams is degraded with increasing x b as K b exp [— (K b x b / F £) m ]. 

The FNS model accounts for the dependence of mechanical response of biological particles 
on particle and indenter geometries 22], Since the information about the particle/indenter 
size is contained in k geometric effects are important in the Hertzian deformation regime. 
When m>l (m<l), the beam failure events become more frequent (less frequent) with in¬ 
creasing x b . For CCMV, 1.8<m<2.1 (Table I), which indicates that the beam-failure rate 
increases with A". When first n beams fail, the compressive force is redistributed among the 
remaining N — n beams, and each surviving beam experiences an increasingly larger tension 
with every next beam failure, which accelerates the failure frequency. The beams fail not 
when f>F£, but under smaller forces F£ ol =F£/ y/em (for 2, F 6 co/ ^0.43F fe *). 

The AFM-based measurements for the empty CCMV shell, empty TrV capsid, full TrV 
virion and full AdV virion are in Figs. S4, S5. Theoretical fit to the experimental average 
FX curves for these particles shows that their deformations are well described by the FNS 
model (Fig. [3]). The obtained Young’s moduli for Hertzian deformation are uniformly smaller 
( 10—100 MPa) than the Young’s moduli for bending deformation (Giga-Pascal range; Ta¬ 
ble I). This correlates with the FX curves being only weakly nonlinear (Fig. [3]). There 
are small variations in the model parameters for the AdV virion due to force application 
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along different symmetry axes. This correlates with onr similar findings for the CCMV shell, 
and implies that the symmetry of local arrangements of capsomer repeats at the point of 
indentation influences its mechanics 22], The failure rates are found to be accelerating with 


deformation (m>l). This fully reflects the conditions of application of compressive force, 
which increases linearly in time speeding up beams’ failure events. Parameters for empty 
and ssRNA-loaded TrV capsids indicate that the difference in particle stiffness is largely due 
to an increase in the Young’s modulus for Hertzian deformation Eh— 0.03 GPa (empty TrV) 
vs 0.14 GPa (full TrV), which suggests that local indentations are resisted in ssRNA-hlled 
particles. These results fit with the previously observed uniaxial deformation of RNA-filled 
TrV into an oblate sphere to maximize the volume available to pack the genome [ll]. Hence, 
confining the large ssRNA genome inside the small particle volume builds internal pressure 
resisting local indentation. 

Previously, the 3D Youngs modulus of the capsid material was estimated using a thin 
shell theory [41, ]l0], [ill, 124]. This assumption is valid for some bacteriophage capsids, but in 
CCMV and TrV the shell thickness is comparable with the virion radius. The FNS model 
accounts for compression of the protein layer under the tip. In the FNS model, the beam¬ 
bending modulus (E b ) is roughly equivalent to the 3D Young’s modulus in the thin shell 
theory. It is estimated at ~0.85 GPa (experiment) and ~0.4—0.5 GPa (simulations) for the 
empty CCMV capsid (Table I). These are similar to yet larger than the values of 0.15—0.30 
GPa obtained with the thin shell theory if, LL0] and 0.28—0.36 GPa from finite-element anal¬ 
ysis (~0.25 GPa) 28], but they disagree with the estimates from several computer modeling 
studies (0.08—0.09 GPa) { 20 , 21], In the study based on spherical harmonics 21], multiple 
deformation modes have also been observed, corresponding to equilibrium deformations of 
the polar regions (tip/surface contact area in FNS model) and the side wall (beams in FNS 
model) of the shell. For the empty TrV capsid, we obtain E b ~ 0.9 GPa (Table I) whereas 
the thin shell theory gives a smaller value of ~0.5 GPa. The lower values of the 3D Young’s 
modulus result from attributing the softer Hertzian deformation mode to bending of the 
capsid shell in the thin shell theory. Indeed, for CCMV and TrV shells, the thin shell theory 
estimates 0.15—0.30 GPa and 0.5 GPa are between the values of AA=0.02—0.03 GPa and 
E b =0.85— 0.95 GPa from FNS modeling (Table I). 

We used parameters of the FNS model and positions of the force maximum at X col =x* H + 
xl from the FX spectra (Figs. [2U3]) to predict the critical force for collapse, F col = f h(x*h) + 
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F™ 1 (xl)—k h^h 1 + K h x * b / Vf/e. The obtained values of F co1 (Table I) agree with their counter¬ 
parts extracted from the FX curves (Figs. [2] [3]), which validates the model. The analytically 
tractable FNS model uniquely combines the elements of continuum mechanics and statistics 
of extremes to accurately describe the mechanical deformation and collapse of biological 
particles. The model can also be extended to characterize the biological particles with other 
regular geometries, e.g. microtubule polymers (cylinder) jlU, 32]. 
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TABLE I. Biomechanical and statistical characteristics determining the deformation and collapse of 
biological particles — CCMV, TrV, and AdV: the Young’s moduli for Hertzian deformation Eh and 
beam-bending E b , the force scale parameter Fj* and shape parameter m. The first (second) entries 
correspond to the exact (approximate) methods of parameter estimation. The model predictions 
for F co1 are compared with the peak forces (in parenthesis) from the spectra (Figs. El [3]). For TrV 
and AdV, the shell thickness was estimated in SM. 


System 


Eh, GPa 

E b , GPa 

F b *, nN 

m 

F co1 , nN 

CCMV (2-fold symmetry; i 

In silico ) 

0.013/0.012 

0.50/0.50 

1.70/1.25 

1.7/1.5 

0.67/0.69 

(0.68) 

CCMV (quasi-2-fold symmetry; in silico ) 

0.011/0.011 

0.37/0.35 

1.50/1.25 

1.4/1.6 

0.58/0.64 

(0.68) 

CCMV (quasi-3-fold symmetry; in silico) 

0.012/0.012 

0.52/0.46 

1.75/1.33 

1.4/1.6 

0.58/0.64 

(0.68) 

empty CCMV (average; in vitro ) 

0.019/0.023 

0.85/0.81 

1.90/1.00 

1.2/1.3 

0.56/0.78 

(0.71) 

empty TrV (average; in 

vitro ) 

0.030/0.036 

0.94/0.81 

1.90/1.1 

1.1/1.2 

0.70/1.02 

(0.69) 

full TrV (average; in vitro ) 

0.140/0.140 

0.95/0.84 

8.10/5.5 

1.1/1.0 

2.91/3.78 

(3.00) 

full AdV (2-fold symmetry; 

in vitro ) 

0.037/0.040 

0.35/0.29 

10.0/5.0 

1.2/1.4 

2.58/4.05 

(3.80) 

full AdV (3-fold symmetry; 

in vitro ) 

0.018/0.019 

0.20/0.18 

11.0/5.0 

1.3/1.7 

3.04/4.15 

(4.30) 

full AdV (5-fold symmetry; 

in vitro ) 

0.021/0.023 

0.14/0.13 

5.10/3.7 

1.1/1.0 

2.03/2.35 

(1.90) 
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FIG. 1. (color) Types of mechanical excitations exemplified using the CCMV shell, (a)-(c) Hertzian 
deformation xh with normal displacements uu p and u par (scheme on (a)) under the influence 
of compressive force (vertical arrow). Dashed contour lines show the tip and particle in their 
undeformed states. Structures in (b) — the native (left) and partially deformed (right) states 
show the amplitude of xh ~3 nm. (c) CCMV shell profile showing parts of the structure with 
high potential energy (> 3 kcal/mol per residue; red) and low potential energy (blue), (d)-(f) 
Beam-bending deformation. The side portion of the structure is partitioned into curved vertical 
beams (top-side view on (d)). Structures in (e) — the partially deformed (left) and pre-collapse 
(right) states reveal the amplitude of 4.3 nm. (f) CCMV shell profile under Hertzian and 
beam-bending deformations showing the potential energy distribution. 
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FIG. 2. (color) Forced deformation of CCMV shell in silico: the average FX spectra (data points), 
obtained from the simulations of nanoindentation (uf=1.0 n m/s, Ru p = 20 nm, and k= 0.05 N/m; 
Fig. S4(a)) along the 2-fold (red), quasi-3-fold (blue), and quasi-2-fold symmetry axes (green), and 
theoretical FX curves (solid lines). The circled bolded X shows the locations of force application. 
We used MD simulations accelerated on GPUs |29(,[30] and nanoindentations in silico of the CCMV 


capsid 


22] (Figs. S1-S3). Structures under the force maxima depict the capsid transitioning from 


the state before the collapse (left) to the collapsed state (right). The top insets: x-based estimation 
of s(X) (left) and dynamics of xh and Xb vs. X (right) in the Hertzian regime I and in the 
transition regime II (dashed and dashed-dotted lines are for the exact method, and solid lines are 
for the approximate piece-wise method of parameter estimation). The bottom inset: schematic for 
piece-wise spectrum modeling; regime I — Xzsxh, and F(X)^Fh', regime II — X=xh + Xb and 
F(X)=F h + F b . 
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FIG. 3. (color) AFM-based forced deformation of empty CCMV and TrV shell, TrV with ssRNA, 
and AdV with dsDNA (top inset which also shows a schematic for piece-wise modeling); see Figs. 
S4(b), S5(a) - (c). For empty CCMV, isj=0.06 /rm/s, Ru p = 20 nm, and k= 0.05 N/m. For empty 
TrV, isf =0.06 fim/s, Rti P = 15 nm, and k= 0.056 N/m. For TrV with ssRNA, isf =0.06 /rm/s, Rti P = 15 
nm, and k= 0.1 N/m. For full AdV with dsDNA, isf= 0.055 fim/ s, Rti p = 15 nm, and k= 0.0524 N/m. 
The average experimental spectra (data points) are compared with theoretical FX (solid) curves. 
The bottom inset shows the 2D surface (Eq.(|5])) constructed with the model parameters from 
CCMV indentation in silico along the 2-fold symmetry axis (Fig. [2j Table I). The red curve is the 
equilibrium average curve F eq X eq with the points formed by the intersection of F(xH, Xb ) with line 
Xb=X—XH (shown for X= 5,7, and 9 nm). The noisy black curves are particular realizations of 
the stochastic FX path. 
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1. Self Organized Polymer (SOP) model of a virus particle 


The SOP model of the polypeptide chain was originally designed to address the mechani¬ 
cal properties of proteins l-4(. The model has been applied to a variety of biological systems 
5-8]. In this work, the SOP model has been used to describe each protein subunit forming 
a virus capsid. In the topology-based SOP model (Fig. IS2lh each amino acid residue is rep¬ 
resented by a single interaction center described by the Ca-atom, and the protein backbone 
is represented by a collection of the C a -C a covalent bonds with the peptide bond length dis¬ 
tance of a = 3.8 A. The potential energy function Usop specified in terms of the coordinates 
of the Ca-atoms {?"*} = rq, rq, ■ ■ ■, rjv (N is the total number of residues) is given by: 


Usop = Ufene + U NB + U J 


' NB 


(SI) 


In Eq. (IS 1|) . the first term is the finite extensible nonlinear elastic (FENE) potential: 


N -1 


tt kR$ ( { r i,i +1 0 , 4 + 1 ) 

Ufene = - 2 ^ — log I 1 -- 


(S2) 


2=1 


where k — 14 N/m is the spring constant, and the tolerance in the change of the covalent 
bond distance is Rq = 2 A. The FENE potential describes the backbone chain connectivity. 
The distance between the next-neighbor residues i and i + 1, is r^+i, and rf i+1 is its value 
in the native structure. To account for the non-covalent (non-bonded) interactions that 
stabilize the native state, we use the Lennard-Jones potential: 
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In Eq. (lS3|) . we assume that if the non-covalently linked residues i and j (\i — j\ > 2) are 
within the cut-off distance of 8 A in the native state, then A^ = 1; A^ = 0 otherwise. 
The value of £h quantifies the strength of the non-bonded interactions. The non-native 
(non-bonded) interactions are treated as repulsive: 
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In Eq. (154|h a constraint is imposed on the bond angle between the residues i, i + 1, and 
i + 2 by including the repulsive potential with parameters £1 = 1 kcal/mol and cq = 3.8 
A. These define the strength and the range of the repulsion. I 11 the SOP model, parameter 
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h sets the energy scale. This parameter is estimated based on the results of all-atom MD 
simulations of the virus particle at equilibrium. 

The dynamics of the virus system is obtained by solving numerically the Langevin equa¬ 
tions of motion for each particle position r; in the over-damped limit: 



dUj(rj) 

dr t 


+ 9i{t) 


(S5) 


In Eq. (1S5D , U t {r^) is the total potential energy, which accounts for the biomolecular inter¬ 
actions ( Usop ) an d interactions of particles with the indenting object — spherical tip ( U tip ; 
see Eq. (lS6l) in section 3 below). Also, in Eq. (fS5l) gi(t ) is the Gaussian distributed zero- 
average random force, and 7] is the friction coefficient. To generate the Brownian dynamics, 
the equations of motion for each C Q -atom are propagated with the time step A t = 0.08r//, 
where t h = CGPl/AbT 1 (A t = 20 ps for CCMV). Here, t l = (ma 2 /£h) 1 ^ 2 = 3 ps, ( = 50.0 
is the dimensionless friction constant for an amino acid residue in water (77 = C m / r L ), 
m~ 3 x 10 22 g is the residue mass, and T is the absolute temperature [J, [dj. To perform 
simulations of nanoindentation of a virus particle, we set T to room temperature and use the 
bulk water viscosity, which corresponds to the friction coefficient r) = 7.0 x 10 5 pN ps/nm. 


2. SOP model parameterization for CCMV shell 


For each virus system, the numerical value of the parameter £h, which describes the 
strength of native interactions (section 1), has been determined from the all-atom MD 
simulations. We performed the all-atom MD simulations of the atomic structural model 
of the CCMV shell (Fig. [5TT) . To obtain accurate parameterization of the SOP model for 
CCMV, we used the crystal structure of the capsid (Viper Data Base; PDB code: 1CWP 
fll | with T = 3 symmetry) in conjunction with the Solvent Accessible Surface Area (SASA) 
model of implicit solvation (CHARMM19 force held) 12], First, we calculated the number of 
native contacts using a standard cut-off distance between the C a -atoms of 8.0 A. The native 
contacts were divided into the inter-chain contacts and the intra-chain contacts. For the 
CCMV shell, there are N intra 20, 554 intra-chain contacts that stabilize the native state of 
the capsid protein, and Ni nter ~ 3,405 inter-chain contacts at the interfaces formed by capsid 
proteins. Next, we calculated the total energy of non-covalent interactions for each contact 
group. The total energy for the intra-chain contacts is E intra = 25, 898 keal/mol and the total 
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energy for the inter-chain contacts is Ei nter = 3, 746 kcal/mol. Finally, to obtain the values 
of parameter Eh we divided the two numbers for each contact group. For the CCMV shell, 
we obtained £i n f rn Ej n tr a l Ty n / ra 1.26 kcal/mol and £ inter Ej r/ j e r /7\/ n f er 1.1 kcal/mol. 
The atomic-level details of biomolccular interactions, contained in these parameters, were 
then exported into the SOP model based reconstruction of the full CCMV particle (Tig. IS2D . 


3. Nanoindentation in silico method 


We employed the methodology of nanoindentation in silico (Fig. 1531) 8], in which the me¬ 
chanical loading of a biological particle is performed computationally. The nanoindentation 
measurements are performed under experimental conditions of dynamic force application 
fit) = rjt (Fig. lS3l) . i.e. in our simulations we use the experimental force-loading rates rj, 
which correspond to the cantilever base velocity Uf = 0.1 — 1.0 /i m/s. Structural transitions 
can be resolved by examining the coordinates of amino acid residues, and biomechanical 
characteristics can be gathered through analysis of the energy output. Because our in sil¬ 
ico “experiment” provides a complete simulation view of particle deformation and collapse, 
it can be used to provide a detailed interpretation and modeling of the experimental force- 
deformation spectra. The full control we have over the system during the nanomanipulations 
in silico can be used to study deformation at different symmetry points on the particle sur¬ 
face, and to relate the force and energy values recorded at any point in the simulation to 
the specific details observed in the particle’s structure. 

In dynamic force measurements in silico, the cantilever base is represented by the virtual 
particle, connected to the spherical bead of radius Ru p , mimicking the cantilever tip (inden- 
ter), by a harmonic spring (Fig. 1531) . The tip interacts with the particles via the repulsive 
Lennard-Jones potential 
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Utip ^ ( - tip 


®tip 
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| Ti Ttip Rtip , 


(S6) 


thereby producing an indentation on the particle’s outer surface. In Eq. (1561) . r, and r tip 
are coordinates of the ith particle and the center of the tip, respectively, £u p = 4.18 kJ/mol, 
and cr t i p = 1.0 Aare parameters of interaction, and the summation is performed over all the 
particles under the tip. For the cantilever tip (sphere in Fig. IS3D . we solve numerically the 
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following Langevin equation of motion: 

dvtip dUtip{vtip) ii o ^ 

V ~df = - dr tip + K " r * p ~ Uft ' ~ nip) ( S7 ) 

where r® ip is the initial position of spherical tip center (yj is the cantilever base velocity; 
k is the cantilever spring constant), and the friction coefficient rj = 7.0 x 10 6 pN ps/nm. 
To generate the dynamics of the biological particle of interest tested mechanically, we solve 
numerically Eqs dSB — (jS5l) for the particle (see section 1) and Eqs. (lS6l) and (lS7l) for the 
indenter (spherical tip). 

The cantilever base moving with constant velocity (uf) (Fig. [S3]) exerts (through the tip) 
the time-dependent force (force-ramp) f (t) = fit) n in the direction n perpendicular to the 
particle surface. The force magnitude, f{t) = r/£, exerted on the particle increases linearly 
in time t with the force-loading rate r/ = KVf. In the simulations of “forward indentation”, 
the cantilever base (and spherical tip) is moving towards the virus capsid. We control the 
piezo (cantilever base) displacement Z, and the cantilever tip position A", which defines 
the indentation depth (deformation). The resisting force of deformation F from the virus 
particle, which corresponds to the experimentally measured force is calculated using the 
energy output. To prevent the capsid from rolling, we can always constrain the bottom 
portion of the particle by fixing select C' a -atoms contacting the substrate surface. 

4. Method of Lagrange multipliers 

To study the dynamical changes in Xh (Hertzian deformation) and xi , (beam-bending 
deformation) and their contribution to the total deformation X = xh + %b (see Fig. 2 in 
the main text), we employed the method of Lagrange multipliers [3|. This method allows 
us to find the values of Xh and Xb that minimize the total deformation force, F(xn,Xb) = 
k H x¥+K b x b s(x b ) (see the inset Fig. 3 in main text), subject to the constraint: X = xn+Xb- 
To that end, we constructed the Lagrange function A(xn,Xb, A) = F(xn,Xb ) + A g(xn,Xb), 
where g(xn, Xb) = Xh + Xb — X and A is the Lagrange multiplier. By calculating the partial 
derivatives of A (xn,Xb, A) with respect to each of the two variables xh and Xb, we obtained 
equations of the form 'V XH , Xb F(xH,Xb) = —\'V XH)Xb g(xH,Xb). Next, by eliminating A we 
arrived at the system of two equations: 

3/2 k H x)( 2 = K b s(x b ) + K b x b s\x b ) (S8) 
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X = x h + x b 


(S9) 


For small deformations x b , we expand the Weibnll survival probability s(x b ) = exp [—{K b x b /F b ) m ) 
in powers of the exponent z = K b x b /F b , and retain the terms up to the first order in z. 
Then, Eq. ( IS8]l becomes: 

3/2 k H x]l 2 - K b ( 1 - (K b x b /F b *) m )(l - m(K b x b /F b *) m ) = 0 (S10) 


Eq. ( IS9|t allows us to eliminate xh by substituting xh 


X — x b into Eq. (IS 10[) above: 


3/2 k H {X - x b ) 1/2 - K b ( 1 - (K b x b /F b *) m )(l - m{K b x b /F *) m ) = 0 (Sll) 


Simplifying Eq. (1S11[) and grouping terms of the same power in x b , we arrive at the following 
polynomial equation: 

aix b m + a-2X b m + a^x b m + a,^x l b + a^x b + clq = 0 (S12) 

which has the following constant coefficients: = rnK b (K b /F b ) 4m , 02 = —2m(l + 

m)K b (K b /F b ) 3m , a 3 = (1 + 4m + m 2 )K%(K b /F b *) 2m , a 4 = -2(1 + m)K 2 {K b /F b *) m , 
a 5 = 9/4 k 2 Hl and a 6 = K 2 — 9/4 k 2 H X. Eq. (S 1S12|) can be solved numerically (for ex¬ 
ample, using Mathematica software) for a given set of parameters kn, K b , F b , and m, and 
for each specified value of the total deformation A". The obtained numerical solution for 
Xh and x b can then be substituted into the expression for F(xn,x b ) (Eq. (10) in the main 
text). 


5. Estimation of the thickness of TrV and AdV with encapsulated genome 

When a virus particle is empty (i.e. it does not contain genomic material) the thickness 
of the beams is equal to the virus shell thickness. The capsid thickness with encapsulated 
genome due to DNA/RNA packing can be estimated assuming that DNA/RNA molecules 
are evenly distributed on the inner capsid surface. The volume occupied by the genome 
considered to form a long cylindrical tube is given by V gen = 7 rR 2 gen L gen , where L gen is 
the genome length and R gen = 1.0 nm and 0.5 nm is the radius of cross-sectional area of 
the dsDNA molecule (for AdV) and ssRNA molecule (for TrV), respectively. The inner 
volume of the virus particle, V)n = V gen = V emp , is the sum of the volume occupied by 
the genome ( V gen ) and the volume of the empty space (V emp ). Here, V in = 4/37 rRf n and 
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Vemp = 4/37r(i£j?i — r g en) 3 , where Ri n = R par — r is the inner radius of the particle shell, r is 
the shell thickness (see Eq. 2 in main text), and r gen is the increase in shell thickness due to 
encapsulated genome. This allows us to express rgen as r gen = R in — (R'f n — ?>/AR 2 gen L gen ) 1 ^. 
For the AdV shell, R in = 36 nm and L gen = 12 //m [lj|; hence, r gen = 2.5 nm. For the TrV 
shell, R in = 11 nm and L gen = 3.06 pm 


15]; hence, r gen = 1.9 nm. These values of r gen for 


AdV and TrV shells were used to estimate Ef, (Table I in main text). 
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FIG. SI. The structure of the Cowpea Chlorotic Mottle Virus (CCMV) (PDB code: 1CWP). The 
side view of the CCMV shell is shown on the right. The protein domains forming pentamers are 
in blue, while the same protein domains in hexamers are in red and orange. The hexamers and 
pentamers, composed of six and five copies of the same protein chain (circled in the black ellipse), 
are displayed on the left. The CCMV capsid is a ~2.8 nm thick shell with a ~26 nm diameter. 
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FIG. S2. Graphical illustration of the coarse-graining procedure involved in construction of a SOP 
model 1, 2] of a polypeptide chain (sections 1 and 2 in SM). Panel (b) shows coarse-graining of the 
atomic structure of the protein subunit forming pentamers and capsomers of the CCMV shell (Fig. 
aSip. Each amino acid residue is represented by a spherical bead of an appropriate radius with 
the coordinates of the Cc-atom (black circles). The protein backbone is replaced by a collection 
of the C a -C a covalent bonds with 3.8 A bond distance. The potential energy function (see Eq. 
(EH)) describes the interactions between amino acids stabilizing the native state of the protein 
chain, and the chain connectivity, elongation due to stretching, and self-avoidance. The coarse- 
graining procedure preserves the secondary structure: a-helices (pink), /3-strands and sheets (blue), 
and random coil and turns (gray). Panel (b) shows the results of coarse-graining of a hexamer. 
Six identical copies of the same protein monomer (coarse-grained in (a)) form a C-based model 
of the hexamer subunit. The hexamers and pentamers are combined to form a coarse-grained 
reconstruction of the full CCMV shell. The SOP model describes the geometry and 3D shape of 
the biological particle. 
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cantilever base 


cantilever tip 



FIG. S3. Schematic of the setup used in nanoindentations in vitro or in silico (see section 3 in SM). 
The biological particle (virus shell) is placed on the substrate. The cantilever base (virtual sphere) 
is moving in the direction perpendicular to the surface of the particle with the constant velocity 
iSf (force-ramp), which creates a compressive force. The force is transmitted to the cantilever tip 
(sphere of radius Ru p ) through the harmonic spring with the spring constant k. The force exerted 
on a particle f(t) = rft (large vertical arrow) ramps up linearly in magnitude with time with the 
force-loading rate rj = KVf, which mechanically loads the particle. The mechanical response of 
the particle can be probed by profiling the deformation force (indentation force) F as a function of 
the cantilever base (piezo-) displacement Z (FZ curve) or as a function of the indentation depth 
X (FX curve). 
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FIG. S4. Nanoindentation of the empty CCMV particle in silico (a) and in vitro (b). Shown 
in different colors for clarity are the FX curves obtained using the cantilever tip velocity Vf = 
0.06 fvm/s (experiment) and Vf = 1.0 /im/s (simulations). In the AFM-based experiments and 
in simulations of nanoindentation of CCMV, we used the cantilever tip with radius Ru p = 20 
nm and the spring constant n = 0.05 N/m. In panel (a), structural snapshots from the left 
to the right, which correspond to the FX curve shown in blue, display the progress of forced 
deformation from the native un-deformed state (leftmost structure), to the partially deformed 
state (middle structures), and finally to the globally collapsed state (rightmost structure). In 
nanoindentation measurements in silico and in vitro , the cantilever tip indents the capsid in the 
direction perpendicular to the capsid outer surface (shown by a large vertical arrow). 
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FIG. S5. AFM-based nanoindentation of the empty TrV particle (a), full TrV particle (encapsulat¬ 
ing the single-stranded RNA genome; (b)) and full AdV particle (encapsulating the DNA genome, 
(c)). Shown in different colors for clarity are the representative force-deformation spectra. The FX 
curves for the empty TrV particle were obtained using the cantilever tip velocity vj = 0.06 /im/s, 
tip radius Ru p = 15 nm, and spring constant k, = 0.056 N/m. The FX curves for the full TrV 
particle were obtained using Uf = 0.06 /im/s, Ru p = 15 nm, and k = 0.1 N/m. The FX curves for 
the full AdV particle were obtained using Uf = 0.055 fin i/s, Ru p = 15 nm, and k = 0.0524 N/m. 
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